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The distribution function of relaxation times in disordered dielectrics has been calculated in the 
random field theory framework. For this purpose, we first consider the dynamics of single two- 
orientable impurity electric dipole in a random electric field E created by the rest of impurities in 
disordered ferroelectric. This dynamics is conveniently described by Langevin equation. Relaxation 
' time r is then a reciprocal probability (calculated on the base of Fokker-Planck equation) of the 

dipole transition through barrier in a double-well potential (corresponding to two possible dipole 
orientations), distorted by a random fields. 

The obtained dependence t(E) made it possible to obtain the expression for relaxation times 
distribution function F(r)(vi& random fields distribution function f(E). Latter function has been 
calculated self-consistently in the random field theory framework. Nonlinear random field contribu- 
tion and effects of spatial correlations between impurities have also been taken into account. It was 
shown that nonlinear contribution of random field gives asymmetric shape of F(t), while in linear 
case it is symmetric. 

Comparison of calculated F(t) curves with those extracted from empirical Cole-Cole (CC), 
Davidson-Cole (DC), Kohlrausch- William- Watts (KWW) and Havriliak-Negami (HN) functions had 
' shown, that they correspond to mixed ferro-glass phase with coexistence of short and long-range 

order. Different forms of F(t) are determined by linear (CC) or nonlinear (DC, KWW, HN) con- 
tributions of random field. 



-a 



> 
o 
in 

(N 

o 



o 
o 



- 1—^ 

X 



I. INTRODUCTION 



The anomalies of dynamic properties are the characteristic feature of the disordered ferroelectrics, polymers and 
composites. In particular, strong frequency dispersion of dielectric or magnetic susceptibility was observed in many 
dipole or spin glasses (see e.g. [1] and ref. therein). The dispersion is commonly attributed to existence in the 
disordered systems of broad spectrum of relaxation times which can be extracted from the observed frequency de- 
pendence of susceptibility [2] . The background of such phenomenological approach is the model of a superposition 
of Debye relaxors with different relaxation times r. In such a model the dynamic quantities like time or frequency 
I ■ dependence of polarization, system susceptibility etc. can be calculated by averaging of single relaxor characteristics 
with relaxation time distribution function F(t) in supposition of parallel (independent) relaxation processes. The key 
point of such approach is the form of distribution function. Only a few simple empirical forms were proposed (rather 
than calculated in some physical model) for F(t). One of them were considered by Frochlich [3] who assumed that in 
arbitrary range [to,ti] F(t) has constant positive value and outside this range it equals zero. But such form of F(t) 
appeared unable to explain the experimental data for real systems with dielectric response described by Cole-Cole 
(CC), Davidson-Cole (DC), Havriliak-Negami (HN) and other more complex forms, than simple Debye equation [4]. 
Different and complex enough empirical forms of F(t) were shown to be necessary to describe non-Debye dielectric 
response [5]. However the extraction of F(t) from the experimental data makes it impossible to understand the 
physical reasons of these properties anomalies. 

Recently the calculations of linear [6,7] and nonlinear [8] dynamic dielectric response were performed for ferroelectric 
relaxors like PMN, PST, PLZT. Random electric field produced by substitutional disorder and other lattice imper- 
fections (which are known to be the characteristic features of all the disordered systems) was considered as the main 
reason of their properties peculiarities via its influence on the barriers between a dipole orientations. Vogel-Fulchcr 
law in relaxation time temperature dependence, stretch-exponential behaviour of dynamic polarization [6] and several 
peculiarities of dc field non-linear susceptibility [8] were obtained. 
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In present work random field theory approach was applied for calculation of relaxation time distribution function. 
The influence of random electric field on the parabolic and rectangular barriers was considered. The relaxation time 
distribution function was calculated with the help of random fields distribution function allowing for linear [9] and 
non-linear random fields contribution as well as spatial correlation effects in the impurity subsystem [10]. 



II. CONNECTION BETWEEN DISTRIBUTION FUNCTIONS OF RELAXATION TIMES AND 

RANDOM FIELDS 

Random field is known to be the characteristic feature of the disordered systems. The substitutional disorder, 
vacancies of ions, impurities and other imperfections are the sources of random electric and elastic field in the system. 
These random fields E make relaxation times r of the dipoles to vary from point to point in a sample, i.e. r = t(E). 
Possible forms of t(E) will be considered below. The distribution of random electric and elastic fields was calculated 
earlier in [9], [10] and [11] respectively in the statistical theory framework. The connection between r and E permits 
to calculate the distribution function of t for known distribution of E by means of theory of probability. In the case 
when t(E) is single- valued (monotonous) function of E the theory of probability gives [12]: 



F(r) = f(E(r)) 



dE(r) 



dr 



(1) 



Here F(t) and f{E) are the relaxation time and random field distribution functions respectively, E(t) is inverse 
function to t(E). The Eq.(nl) can be thought about as simple change of variables in normalization integrals 



J F{r)dr = J f(E)dE 



(2) 



where integration is performed over domains of r and E variation. The only difference is the modulus of the derivative 
because the distribution function has to be positive. 

In more general case of non- monotonous behaviour of t(E), when one r value correspond to several E- values 
(Ei, E%...E n ), the space of E values should be divided into n regions (containing E\, Ei-..E n points), where function 
t(E) is monotonous (see Fig.l). Now, for the entire .E-domain, F(t) can be represented as a sum of expressions like 
Eq.(pl) over regions Ei(r) of monotonous behaviour of t(E) [12]: 



F(r)=J2m(T)) 



dEi(r) 



dr 



(3) 



Eq. ([|) is general expression for the distribution function of one quantity via the distribution function of another one 
for a given connection between them. In particular, Eq. (|^) can be used for the calculation of random field distribution 
function f(E) allowing for nonlinear and spatial correlation effects via that in the linear case. In general case internal 
random field E(r{) can be written as: 



E~f(rt) = e 7 (rt) + 



E 



l (rt) 



(4) 



£ 7(^) = ^2(£k)y(rij) (5) 

where (£fc) 7 is 7-component of the field produced in the observation point r| by a source of k th type (e.g. dipoles, 
point charges, dilatational centers) situated at the point rj; a m are the coefficients of nonlinearity of m th order of 
the host lattice, their dimension being the inverse electric field in (to — l)-th power. 

In the linear case (a m — 0) the random field distribution function /; can be calculated in the statistical theory 
framework, namely 

fi(e) = 6(e-e'( rij )), (6) 

where averaging (marked by bar) of (5-function makes it possible to count all the spatial configurations of random field 
sources e'(rij) leading to definite random field value e 7 = e. Averaging, neglecting the correlations between random 
field sources in Eq.(||), leads to well-known expression of statistical theory [13,14]: 
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n k F k {t))dt 



(7) 



F k (t) 



tfV(l-exp(ite fe (r))) 



(8) 



Here n k and £&(r) are respectively the concentration and electric field of k th random field source. 

One can see, that the form of F k (t) and thus fiie) depends strongly on the form of £fc(r). So, fi(e) could not only 
be of Gaussian form (which is realized, e.g., for point charges with e k (r) ~ \/r 2 ), but also of Lorentzian form (this is 
the case for electric dipoles with e k (r) ~ 1/r 3 ). 

In the disordered dielectrics the main source of random fields is impurity electric dipole which can have several 
discrete orientations in a host dielectric lattice (see Ref. [9] for details). In this case, which is of interest for 
present consideration, the function fi(s) chould be calculated with additional quantum statistical averaging over 
aforementioned discrete orientations. This was done self-consistently in Ref. [9] (see also [10]) for particular case of 
two-orientable dipoles, point charges and dilatational centers as the random field sources: 



fl(E) 



1 

2^ 



exp 



it(E - EqL) - mBi \t\ - n 2 B 2 \t 



3/2 



n 3 B 3 r 



dt 



(9) 
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L = / fi(E)t&nh((3E)dE. 



Here L = ((d*))/d* is long-range order parameter (the number of coherently oriented impurity electric dipoles, see 
Ref. [9] for dtails), Eq — A-K(n 3 d* 2 ) / Eo is the mean value of random electric field (in energy units), produced by the 
impurity dipoles, d* — 1/Sdj (sq — 1) is the electric dipole moment, 7 and eo are reaspectively the Lorentz factor 
and static dielectric permittivity of the host lattice, ni, n% and n 3 are respectively the concentrations of dilatational 
centers, point charges and electric dipoles, Ze and Oo are the point defect charge and elastic moment, p and v are 
the host lattice piezoelectric component and Poisson coefficient respectively. 

Thus Eqs. (6), (7), (8) give the possibility to calculate the distribution function of random field induced indepen- 
dently, i.e. without correlation, by any number of random field sources. 

Spatial correlations and nonlinear effects can be taken into account on the base of Eq. (3) (a n 7^ 0) . For such more 
general case the random field distribution function can be expressed through /;(e) with the help of expression like 
Eq.(2) [15]: 



de n {E) 



dE 



where e n are the real roots of the algebraic equation 

E - e - a 2 e 2 







(10) 



(11) 



ds n (E) 
dE 



We have to emphasize that even for simple enough (e.g. Gaussian, Lorentzian) form of fi(s), coefficient 

in Eq.(9) will be function of E by virtue of Eq (|TT|) . So, the nonlinear distribution function f(E) will never be the 
simple algebraic sum of fi{e n ). The shape of f(E) was calculated recently with the help of Eqs. ([10|), ( |TT| ) for the 
contribution of nonlinear terms of the second (a 2 ^ 0) and third (0:3 ^ 0) order [15,16]. 



III. INFLUENCE OF ELECTRIC FIELD ON RELAXATION TIME 



The ordinary ferroelectric materials like PbTiOa, BaTiOa with polar long-range order are known to have unique 
relaxation time characterizing the rate of macroscopic polarization restoration after external perturbation of the 
system. 
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Disordered systems like relaxors have broad spectrum of relaxation times, because their most probable states are 
dipole glass (DG) with short-range order polar clusters of r c size, imbeded into paraelectric phase and/or mixed ferro- 
glass phase (FG) with coexistence of short- and long-range order. Since an orientable electric dipole is usually the 
nucleation centre of a cluster, the random reorientation of the dipoles can be the main mechanism of the relaxation 
in the aforementioned disordered systems. It is commonly belived that the probability of a dipole reorientation 
temperature dependence is obeyed to Arrhcnius law 

r = - eXP {-T) (12) 

where U is a height of a barrier between equivalent dipole orientations. 

To obtain the connection between relaxation time r and random field E 7 let us consider the relaxational dynamics 
of single impurity (for definiteness two-orientable) dipole in a random electric field, created by the rest of impurity 
dipoles. It is clear, that in such relaxation process the random field will be time dependent and will play a role of a 
random force, exerted on the given impurity dipole. For such motion we can write following Langevin equation (see, 
e.g. 0) 

where x is the displacement of impurity ion (so that impurity dipole moment d = ex, e is ion charge), U (x) determines 
the actual potential energy of dipole (e.g. for two-orientable dipole U(x) is a conventional double-well potential) and 
f(t) is a random force with known correlator. Without loss of generality we can put 

< f(t) >= 0, < /(*)/(*') >= DS(t - 0, (14) 

where D is diffusion coefficient (see below). Langevin equation ([l3]) is nonlinear stochastic differential equation. It 
can be used, in particular, to determine the probability for impurity dipole to transit from shallower to deeper well of 
double- well potential (i.e. overbarrier "bounce"), deformed by a random field (Fig. 2a). This problem is equivalent 
to the problem of Brownian particles diffusion through barrier in the potential U(x) (see [l8| and refs therein). For 
the probability n(x, t) that diffusing particle has coordinate x at time t the Fokker-Planck equation can be derived 
from (|l^). For two-orientable dipole (one-dimensional case) it reads 

dn dU(x) 
j = -D- on— . 

ox ox 

Here D is diffusion coefficient for particles from well A to well B of potential U(x) (Fig. 2a). For it the Einstein 
relation holds 

D = bT, 

where b is particle mobility, j is diffusion flux. It is seen that dn/dt — in Eq. (|l^) when 

n = n eq (x) = Acxp (^-^jr^J , 0-6) 

i.e. when n is equilibrium Maxwell-Boltzmann distribution. 

It can be shown (see, e.g. and refs therein) that probability of transition from well A to well B (i.e. reciprocal 
relaxation time) is ratio of stationary (at dn/dt = 0) flux jo and number of particles Na in well A 

W < 17) 

Calculation with respect to boundary condition n(B) — (reflecting the fact that at t = all particles are in well A) 
yields 

j = - , N A = /.(A)exp (m) dx (18) 

A 
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One can see from Eqs. jl7|), (18), that the relaxation time depends mainly on the form of potential U{x). It 
can be shown that the main contribution to the integrals originates from potential U(x) maximum region. So, we 
can approximate U(x) by parabola near this point and integrate from — oo to +00 (this is standard steepest descent 
method for integrals calculation). For zero random field (E = 0) U(x) = Uq(x) is symmetric with respect to y axis 
and this procedure gives standard Arrhenius law (^) for relaxation time. For nonzero random local (i.e. in the given 
impurity location) electric field the symmetric form of Uq(x) is distorted 

a 2 

U(x) = Uo(x) ± eEi oc x, U = U-ax 2 + bx 4 , U = —. (19) 

Here U is barrier height and for undistorted potential we choose the simplest polynomial approximation. In subsequent 
calculations we will call this barrier shape "parabolic". 

Consider now aforementioned potential Uq(x) = U — ax 2 + bx 4 , which has maximum at x\ — and minima at 
X2 : 3 = ±y / a/2b. In the spirit of steepest descent method we considered the approximation of U{x) ( |l^ ) near its 
maximum at x = for small fields E. This made it possible to carry out the integration in Eqs. (15), (16) by the 
steepest descent method, that yields 

P= l = L exp (- U +^ 2 /* a ) (20) 

where 1/tq is temperature and field independent combination involving constant b and Uq{x) parameters. Note, that 
at E = we once more obtain Arrhenius law from (p0|). 
For rectangular potential of the form (see Fig. 2b) 

0, — xo — A < x < —xq 

Unlr) = { U, -X <X<X (21) 

0, xo < x < xo + A 

the integrals (0), ( |l8| ) can be calculated exactly, which yields 

p = a 2 Dey^{~axo) 

[2/aexp(U/T) sinh(aa;o) + 4/a sinh(aA/2) cosh(axo + aA/2)] (exp(aA) - 1) ' 1 ' 

a = — , d = ex - 

It is seen from the Fig. 2a that actually A -C xo- In this case we can substantially simplify Eq ( p2| ) and obtain 

1 1 / U±dE loc \ 

- = -exp^ _j (23) 

It is seen that at Ei oc — we again arrive at Arrhenius law. Note that steepest descent method gives also Eq.(^). 



One can see from Eqs. (|20|), (23), that electric field decreases the barrier for dipole orientation along field direction 
and increases it for opposite one. Thus, random electric field increases relaxation time. 

The final form of relaxation time dependence on E can be obtained after averaging Eq.(21) over possible orientations 
of electric dipole d. Quantum-statistical averaging with the energy in the form 3? = —dEi oc gives 

spexp(d*S/r-3?/T) ch(2d*E/kT) fU\ 



In Eq. (24) we substituted product dEi oc by d*E where effective dipole moment d* = (sq — l)'yd/3 (eo and 7 are 
respectively host lattice permittivity and Lorentz factor). 

Note, that such type of electric field influence on the barrier height was supposed by us earlier (see e.g. [6-8] and 
ref. therein). 



In what follows we will perform the calculations for rectangular barriers on the base Eq.(24), because the influence 
of electric field on parabolic barrier (n9) is qualitively the same. 
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IV. RELAXATION TIME DISTRIBUTION FUNCTION 



A. General equations 

Eq.(2) with respect to Eq.(22) makes it possible to calculate relaxation time distribution function. In particular 
Eq.(22) leads to the following connection between electric field and relaxation time 



E(t) = fcTarccos/i (~(t + \/t 2 + 8) J = kTf„±(t);t = =;t>l (25) 

where sighs " ±" correspond to two branches of the function 

arccos h(x) = ln(x ± \/ x 2 — I) (26) 

Allowin g fo r these two branches as i = 1,2 in Eq.(2) and substituting E(t) for E in Eq. (^) (linear case) or in 
Eqs.([lO|), ( jll| ) (nonlinear case) one can obtain general expressions for relaxation times distribution function. Since 
these expressions appeared to be quite cumbersome, we considered the case when only electric dipoles are the sources 
of random electric field i.e. n\ — n 2 — in Eq.(§). In this case Eq.([}J) yields 

ME) = ^exp (- (E ~^ L)2 ) , A = v^ft, (27) 

i.e. linear random field distribution function fi(E) has Gaussian form. Eq.(25) with respect to Eqs (||) and ( p5| ) leads 
to the following expression of relaxation time distribution function 

exp[-b 2 {vfo-(t)-L(v,z)f]X (28a) 



\/ir 



dt 



exp \-b 2 (vf 0+ (t) -L(v,z)) 2 ] + 



dfp- 
dt 



dfo± =± _1 



dt \/2 y/t 2 + tVW+8 - 4 

where we introduced the dimcnsionlcss parameters suitable for numerical calculations 

E 



h = 



kT 



2%AOT 



= vbo ; bo = 



1 



= — VTBttz: 

2V^Bl 2 



(28b) 



(29) 



kT 
~E~ 



z = n 3 r c ;v 

Let us check the normalization of F(t). The normalization condition for f{E) 



oo 

j f(E)dE = l 



(30) 



transforms for F(t) into contour integral, where contour runs first from oo to 1 over t axis and then back - from 1 to 
oo (see Fig. 3). 

To perform integration we should pass from integration over dt to the integration over dfo (with respect to /o±(l) = 
0, /o+(oo) = oo, /o-(oo) = -oo). This yields 



exp |-(6i/o - b 2 L(v, z)) \dfo+ / exp I — (6i/q - b 2 L(v, z)) \ df 



(31) 



exp 



J -b 2 L(v,z)) 2 ]df = l 



Therefore Eqs.(|28a|) 
ing function F(t) may a 



represent normalized distribution function of relaxation time. Any other integral involv- 
so be evaluated by this procedure. 
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One can see that F(t) depends on random field sources concentration and temperature via dimensionless parameters 
z, v and order parameter L(z, v). The dependence L(z, v) actually determines the phase diagram of the system under 
consideration. We calculated it earlier both for linear [9] and nonlinear [10] random field contributions. It was shown, 
that L = for 713 < n 3ci T < T g , where n 3c and T g are respectively critical concentration and freezing temperature 
for dipole glass (DG) state; L — 1 for z — > 00 (mean field approximation) and T <C Eq = kT cm f, where T cm f is 
ferroelectric (FE) phase transition temperature. In intermediate range of concentrations (00 > n 3 > n 3c ) only the 
part of dipoles is ordered, i.e. < L < 1 and ferro-glass (FG) phase with coexistence of long and short range order 
appears. 



Thus, Eqs.(E8a|), (28b) permit to calculate the distribution function of relaxation times for DG, FG, FE and 



paraelectric phases via L(z, v) dependence in these phases. 

B. Relaxation time distribution function in mean field approximation 

In the ferroelectric and paraelectric phases all properties of disordered dielectrics, including distribution of relaxation 
times, can be calculated in mean field approximation. In this approximation random field distribution function has 
the form of 5-function, i.e. f m f(E) = 8(E — EqL) (linear case) and 6(E — EqL(1 + a^E^L 2 )) (nonlinear case [10]) 
with L 7^ for FE phase and L = for paraelectric phase. Substitution of these functions into Eq.(l) gives e.g. for 
linear case 



F mf (t) = 5(E(t) - E L) 



dE{t) 



df 



(32) 



Taking into account that f(x)S(x — a) = f(a)5(x — a) and 8(f(x)) — 3( x — x k)/ \f'{ x k)\ (see, e.g. ,[18]), where Xk 

k 

are the real roots of equation f(xk) = 0, one can rewrite Eq.(p2) in the form 

F mf (t)=S(t-t mf ) (33) 

where t m f was shown to be the unique root of the equation E(t) — E L = (linear case) or E{t) — E L(l+a 3 EQL 2 ) = 
(nonlinear case), see Fig. 3. One can see from Eq. ( p5| ) that in paraelectric phase (L — 0) t m f = 1, i.e. r m f = ni = 
tq exp(U/kT), whereas in FE phase (L ^ 0) additional temperature dependence of relaxation time appears. In 
particular for LE /kT > 1 r m / = t exp((f7 + LE )/kT) (linear case) or r m f = t exp[([/ + LE (1 + a 3 E 2 L 2 ))/kT] 
(nonlinear case), because of mean field influence on barrier height. 



C. Numerical calculations of relaxation time distribution function 



The calculation of relaxati on t ime distribution function beyond of mean field approximation were performed nu- 
merically on the base of Eq.(28a) (linear case) and on the base of Eqs. (|J), (|ic|), (fill), (HI) (nonlinear case). Since 
the spectrum of relaxation times in disordered systems is very broad (1 < t/tq =t< 00), we, following literature on 
this subject, plotted F(t) in logarithnic scale in t. By the same reasom we did not include dfo±/dt ~ d(lnt)/dt into 
our plots. Also, fi(E) was always taken in Gaussian form (p7|). 

The results of calculations of relaxation time distribution function in the linear case are represented in Fig. 4. To 
exclude a paraelectric phase we choose v = T/T cm t < 1 and considered two values of dipole concentrations z = nr~ = 1 
and 10 (curves 1 and 2 respectively). Due to the values of parameters chosen, curve 1 corresponds to F(t) for FG 
phase whereas curve 2 illustrates the transformation of F(t) into the form close to 5-function, which is realized as 
z increases (system approaches to FE phase). The shape of curve 1 reflects the Gaussian form ( p7^ of fi(E) . For 
non-Gaussian shape of fi{E) (e.g. Lorentzian, Holtzmarkian etc, corresponding to different random field sources), 
the shape of F(t) in linear case will also reflect them. Since all aforementioned curves are known to be symmetric 
functions, the linear relaxation time distribution function also must be symmetric. 

Let us proceed to calculation of F(t) in nonlinear case. Since the solutions of Eq. ( pd| ) depend on the lattice 
symmetry we considered a host lattice with a centre of inversion in the paraelectric phase, that is characteristic for 
many disordered systems. In such a case only the terms with odd powers of e remain in Eq.([ll]). In our calculations we 
retained only the first nonlinear term in Eq.([ll]), i.e. o^e 3 . The sign and value of nonlinearity coefficient a 3 = aoE^ 
appeared to influence strongly the shape of F(t) (see Fig. 5 for ao > and Fig. 6 for ao < 0). One can see that for 
ao > F(t) broadens and shifts towards larger t with ao increase. Such behaviour can be the result of barrier increase 
due to nonlinear random field contribution. Shape of F(t) transforms from almost symmetric to slightly asymmetric 
with ao increase. For example, for ao = 1 the right "wing" of the curve is just Gaussian while left one decays faster 
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then Gaussian. Strongly asymmetric shape of F(t) is peculiar feature for ao < 0, the line asymmetry being larger at 
larger |cko | (see Fig. 6). The narrowing of F(t) and its maximum position shift to the smaller t with |ao| increase can 
be the result of the barrier decrease. More essential asymmetry of F(t) for ao < in comparison with ao > is the 
result of presence of only one root of Eq. ( |TT| ) (and symmetric shape of f(E)) for ao > 0, whereas for ao < there 
are three different real roots of Eq. JTl|), which breaks f(E) symmetric shape (see [15]). The curves depicted in Figs. 
5 and 6 correspond to FG phase with coexistence of short- and long-range order. Calculations of F(t) for DG state 
(where only short-range ordered clusters exist ) lead to very broad curve with very slow decay even in ln(t) scale. 



V. DISCUSSION 



The distribution of relaxation times is usually extracted from observed frequency dependence of dielectric suscepti- 
bility. This dependence for ordered systems is known to be described by Debye law with the unique relaxation time. 
To describe the observed dynamic susceptibility of the disordered ferroelectrics, polymers and composites, several 
empirical functions were proposed as of Debye law generalizations (see e.g. [4]). Among them the most known are 
following: 



1 + (iuJTcc) 1 - 11 ) (a) 
(I + ilutdcV 13 (b) (34) 

11 x 1 (l + (iur HN y)- S (c) 



Eq.@ (a), (b), (c) represent, respectively, Cole-Cole (CC) (0 < k < 1), Davidson-Cole (DC) (0 < < 1) and 
Havriliak-Negami (HN) (7 < 1, 5 < 1) functions. 

All these functions are written in the frequency domain, whereas in time domain Debye relaxation is usually 
generalizes to Kohlrausch- Williams- Watts (KWW) relaxation function, which is also currently called a "stretched 
exponential" 

$(t) = exp(-t/ Tww ) a ,0 < a < 1 (35) 

Since there is no analytical expression for Eq. ( |35j ) in frequency domain, the numerical calculation of its Fourier 
transform was performed [4]. It was shown that at 7<5 = a 1,2 ' 3 the results are identical to those for H-N function. 

Empirical laws ( |34]) , ( |35| ) were applied for many years for the description of slow relaxation processes in conventional 
glasses, polymers, composites, disordered ferroelectrics etc. The data obtained by several experimental techniques, 
including dielectric spectroscopy, nuclear magnetic resonance, quasielastic neutron scattering , kinetic reactions etc. 
were sucsessfully fitted by these laws (see e.g. [4], [19]). It was supposed , that physical origin of the laws ([54]) and 
( |35| ) was some distribution of relaxation times F(t), so that 

00 

£ * M ~ £ °° = / T- —F(r)d(lnr) (36a) 

£0 - £00 J1+ 



00 

exp (— t/rww) a = J exp (— t/r W w) F(r)d(lnr) (36b) 


The expressions ([34]), ( |36a| ),(??) made it possible to extract relaxation time distribution function for all the empirical 
laws. The results of such extraction, obtained in Ref. [5] are represented in Fig. 7 for Debye CC (x = 0,2), DC 
(/3 = 0, 6), KWW (a = 0, 42). Due to aforementioned relation between KWW and HN laws the shape of F(t) for HN 
is similar to that for KWW function (see Fig. 7). The obtained distribution functions are empirical ones with strongly 
different shapes. Namely, F(lnr) can be symmetric (CC), asymmetric (KWW, HN) and strongly asymmetric (DC) 
function. To the best of our knowledge, the physical mechanisms of such behaviour is not known out up to now. 

To clarify the physical mechanisms let us compare the empirical curves from Fig. 7 with the curves in Figs. 4, 5, 6, 
obtained on the base of random field distribution function. As it was shown in Section 4 the symmetric form of F(r) 
is pertinent to linear approximation what is valid for small enough random field sources concentrations. With increase 
of random field sources concentration nonlinear and correlation effects become essential. They result in F (In t/tq) 
small asymmetry for ao > (Fig. 5) and strong asymmetry for ao < (Fig. 6). The comparison of the curves in Figs. 
6 and 7 shows that the shape of the curves 2 and 3 in Fig. 6 looks like that for KWW and DC laws respectively (Fig. 7). 
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Relative position of their maxima is in conformity with that for curves 2 and 3 in Fig. 6. The coinsidence of tcc and 
Debye relaxation time td (see Fig. 7) can be expected also from calculated curves in Fig. 4 because with increase of 
random electric dipoles concentration curve 2 has to be transformed into 5-function (see Section 4). Since calculated 
curves in Figs. 4, 5, 6, correspond to FG phase with coexistence of short- and long-range order, we can conclude that 
this coexistence as well as small (CC law) or large nonlinear effect contribution with negative coefficient (KWW, DC 
laws) are the physical background of the considered empirical laws. Note, that nonlinear and correlation effects with 
negative nonlincarity coefficient decrease the system order, i.e. they makes it "more disordered", as it was shown 
recently [10]. Under such conditions the role of random fields and thus the distribution of relaxation times becomes 
more essential in the peculiarities of physical properties of disordered dielectrics. 

VI. CONCLUSION 

The coincidence between the relaxation time distribution function obtained from empirical laws and calculated on 
the base of random fields distribution function gives evidence that random fields are indeed the main physical reason 
for the distribution of relaxation times. Variety of distribution function forms is due to "degree of disorder" (i.e. 
presence of different amounts of different random field sources) in the systems under consideration. 

Suggested general formalism for relaxation times distribution function calculation has pretty general nature (because 
of generality of statistical method, see |l3|,[l4]]) within assumptions made. It can be easily generalized to other 
disordered systems like spin glasses. 
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FIG. 1. Schematic plot of nonmonotonous behaviour of t(E) function. One value to corresponds to four E values (Ei — E4). 
In the regions 1-4 function t(E) is monotonous (it decreases in regions 1 and 3 and increases in 2 and 4). 

FIG. 2. (a) - Parabolic barrier in the external electric field E. It is seen that barrier height at E ^ is larger than that at 
E = 0. Parabolic approximation for steepest descent method is also shown, 
(b) - Rectangular barrier at E = (solid line), E < (dashed line) and E > (dotted line). 







FIG. 3. Schematic plot of dependence E(t) (|25| ). Arrows show the integration contour in (|3l[). We also show the unique 
root t m f of equation E(t) — E m f = (see Eqs (|3*4|, (|35[)), i5 m / = EqL (linear case) or E m f = E L(l + aoL 2 ). For the sake of 
illustration E m f is plotted only for the case L > and ao > 0. 



FIG. 4. Relaxation time distribution function for small random field contribution - linear case. Parameters of calculations: 
v = 0, 3, z = 1 (curve 1), z = 10 (curve 2). 

FIG. 5. Relaxation time distribution function for large random field contribution - nonlinear case, with positive (a) and 
negative (b) coefficient of nonlinearity. Parameters of calculations: v = 0, 5; z = 2; ao = 1; 0, 5; 0, 1 (curves 1, 2, 3 respectively 
in (a)); v = 0, 3; z = 1; ao = —0, 01; —0, 1; —0, 3 (curves 1, 2, 3 respectively in (b)). 



FIG. 6. Relaxation time distribution function for Debye (D), Cole-Cole (CC), Davidson-Cole (DC) and 
Kohlrausch- Williams- Watts (WW) laws [5]. 
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